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Extending recent work of Corrado, we derive an algorithm 
that computes rigorous upper and lower bounds for rectangle 
scan probabilities for Markov increments. We experimentally 
examine the closeness of the bounds computed by the algorithm 
and we examine the range of tractable input variables. 



1. Introduction. Let n balls randomly fall into d boxes, each ball with 
probability pi into box i G {1, . . . , d}, independently from all the other balls. 
What is the probability that there exist i adjacent boxes in which together 
lie more than k balls? Formally, if we turn to compute the probability of the 
complement: Let N ~ M„p be a multinomially distributed random variable. 
Task: Compute 

F{N^ + ... + Ne<k,..., Nd-e+i + ■ ■ ■ + < k) 

In this paper, we derive an algorithm that allows fast computation of this 
probability. 

Such probabilities are needed as p-values for tests that check data on clus- 
ters. For example: Let n = 500 patients arrive at a clinic in d = 365 days. 
We compute the probability that there exist three successive days in which 
together more than 15 patients arrive. From the line for = 15 in Table Oon 
page fm below, we get the approximate value 1 — 0.9979961 = 0.0020039 with 
an absolute error less than 10"''. As this probability is so small we would, 
if the described event occurs, reject the hypothesis that the patients arrived 
independently and hence suspect that there must be a reason for this cluster. 

The support D = {x G Nq : Xi + . . . + Xd = n} of the multinomial 
distribution M„^p is finite. Hence we could compute the desired probability 
as follows: For each x E D with xi + . . . + xe < k, . . . , Xd-e+i + . ■ ■ + Xd < k 
compute the probability P(A^ = x) = n!/(xi! . . .Xd^-)Pi^ ■ ■ -P^ and sum up 
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these values. But because the support D is large, this procedure takes much 
time. For example: If n = 15, c/ = 12 it took 4 seconds to compute the 
probability vr := P(iVi + N2 + < . . . , Nd-2 + A^d-i + A/'d < 5) on a 3 GHz 
desktop pc, for n = 15,d = 25 it already took 8-9 hours. 

To derive a faster method in Sections 2 and 3 of this paper, we use a fact 
already utilized by Corrado [3], namely that the multinomially distributed 
random variable is a Markov increment, see Section 4. In this paper, a 
Markov increment is a vector {Yi, . . . ,Yd) of discretely distributed ran- 
dom variables with values in a group {X, ■) with the property that {Yi,Yi ■ 
Y2, . . . ,Yi ■ ■ ■ Yd) is a Markov chain. Our method actually works for Markov 
increments in this generality. For example, the computation of the probability 
vr for n = 15, d = 25 with the new method takes less than one second. 

In Sections 5 and 6 we turn to computer-implementations of our algorithm 
within the IEEE-754-standard [2] for floating point computer arithmetic. 
The floating point number systems according to the IEEE-754-standard that 
usual computers work with have the following properties: The exact result of 
an operation on two floating point numbers, e.g. addition, need not be a float- 
ing point number again. In that case, the computer returns a floating point 
number that is as close as possible to the exact result. The difference between 
the returned value and the exact result is called rounding error. Because 
of rounding errors, computed values, e.g. probabilities, are usually just ap- 
proximations for the exact values and the goodness of the approximation is 
not known. One can switch the rounding mode of the machine in such a way, 
that in every operation it returns the minimal floating point number which is 
greater or equal than the exact result. This "rounding up" mode can be used 
to compute upper bounds for the exact value, if only positive numbers occur 
and only additions and multiplications are performed. In the same way, per 
"rounding down" mode, lower bounds can be computed. Thus one gets an in- 
terval whose bounds are floating point numbers and in which the exact value 
is known to lie. The accuracy of the approximations can easily be estimated, 
because the two bounds of the interval are known. In Section [6l we present 
an implementation of our algorithm within R. For definiteness, we assume 
that all computations are done in double-precision according to the lEEE- 
754-standard. We analyze the accuracy of the R-implementation and compare 
it to the best possible accuracy in lEEE-Double-Precision-computations of 
probabilities, which we examine in Section 5. 

To sum up: This work extends Corrado 's by clarifying the underlying 
Markov increment structure, by allowing the computation of scan proba- 
bilities and by providing rigorous numerical bounds. 
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2. An algorithm that computes rectangle probabilities for Markov 
increments. We derive an algorithm that computes rectangle probabilities 
for Markov increments. It is based on the following recursion formula: 

Theorem 2.1. Let Y — (1^)^=1 be Markov increment of a Markov chain 
{Xk)'l^i which takes values in a group {X, •) . Let Ai, . . . , G X be countable 
sets. Then the probabilities 

p{k, x) := F{Xk = X, Yi e Ai, . . . , n e Ak) 
for k & {1, . . . ,d} and x & X fulfill the recursion 
(1) p{k, x)=Y^ ^i^k - X I Xk-i = xy-')p{k - 1, xy-') 

for k > 2. Here and throughout, we use the convention ¥{A\B) — ¥{A fl 
B)/F{B) :=0 z/P(S) =0. 

Proof. The functions f^ : X"^ —f X defined by fk{xi,X2) = x^^X2 have 
the property that Yj. = fk{Xk_i, X^) and fk{-, x) is bijective for every x G X. 
Using this (which is actually all we need, so the method works not only for 
Markov increments but actually for any functions of two successive states of 
a Markov chain having the above bijectivity property) and writing gk{x, •) := 
fk{-,x)-^, we get: 

F{Xk^x,YieA^,...,YkeAk) 
= J2 = x,Yk = y,Yie Ai,..., n_i e Ak-i) 

= J2 H^k = X, Xk-i = gk{x, y),Yie Ai,..., Yk-i e Ak-i) 
yeAk 

= J2^i^k^x\Xk-i^gk{x,y)) 

y&Ak 

X P(Xfe_i = gk{x, y), Fl e ^1, . . . , n_i e Ak-i) 
In the last step the Markov property was used. □ 

Prom the recursion formula we can derive the following algorithm that 
computes the probability P(Yi G Ai, . . . ,Yd e Ad). Let Ai, . . . , be finite, 
so that we get a finite algorithm. 
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Algorithm A: 

1. For every x & Ai compute the value p{l, x) = P(Xi = x) 

2. For every k G {2, . . . , d}: 

For every x G Ai ■ . . . ■ Ak compute the value p{k, x) with formula ([T]) 

3. Compute 

xeAi-...-Aa 



Here, let Ai ■ ■ ■ An := {ai ■ ■ ■ a„ : ai G Ai, . . . ,an G An}, if A" is a group and 

Ai,...,AnCX. 

3. Computing rectangle scan probabilities for Markov increments. 

In this section we describe how to compute a rectangle scan probability 

q := P{Yi ■ ...-YeeA^,..., Y^.i+i ■ . . . ■ e A^-e+i) 

for a Markov increment Y. 

We use the following obvious and well-known lemma: 

Lemma 3.1. Let X be a countable set and {Xk)f^^ an X -valued Markov 
chain. Let Wk ■= (X^, . . . ,Xk+i-i). Then {Wk)'lz^{^^ is an -valued Markov 
chain with transition probabilities 

¥{Wk+i =w\Wk = v) = = I Xk+i-i = V,) 

for v,w E X^ with f{Wk = v) > and V2 = Wi, . . . ,Vi = 

The desired rectangle scan probability for the Markov increment Y can 
be written as a rectangle probability for the increment V of W: If we set 
Bk ■= {{yi, ■ ■ ■ , ye) e X^\yi ■ ...-ine At} we have 

q = F{VieBi,...,Vd^e+ieBa-e+i) 

because Vk = {Xk - Xk-i, . . . , Xk+e^i - Xk+e) for /c G {2, . . . , - £ + 1}. 

The sets Bi, . . . , Bd~i+i ^'^^ possibly infinite so the Algorithm A from the 
last section would not work. But if there exist finite sets Mi, . . . , Md^e^i C X^ 
with 

P(\/i G 5i, . . . , Vd-i+i G Bd-i+i) = G Ml, . . . , Vd-e+i G Md-e+i) 

we can apply the Algorithm A and thus are able to compute the desired 
probability. 
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Example: U X = (Z, +) and F is a Markov increment with Yi, . . . ,Yd > 0, 
then for finite sets Ai, . . . , Ad-^i+i C Z the probability 

P(ri + . . . + e Ai, . . .,Ya-i+i + ... + YaeAa) 

equals 

P((ri, . . . , F,) e Ml, . . . , . . . , y,) e Ma-,+i) 

with Mk := {(z/i, . . . , ?/^) e N^lyi + . . . + e A^}, which are finite. 

4. Examples for Markov increments: Multinomially and multi- 
variate hypergeometrically distributed random vectors. By hn,p{k) — 
(^)p'^(l — pY'^ we denote the binomial density with parameters n e N and 

p e [0, 1]. By h.n,r,b{k) = Q in-k) / {''ti^) denote the hypergeometrical den- 
sity with parameters r,b eNq and n E {1. . . . .r + b}. 

Multinomially distributed random vectors as well as multivariate hyper- 
geometrically distributed random vectors are Markov increments, hence the 
results from the last two sections are applicable in these cases. More precisely, 
we have the following two propositions, as easy calculation with density for- 
mulas and cancelling yield. 

Example 4.1. Let {Ni,...,Nd) ~ M„^p be a multinomially distributed 
random variable and '■— X]j=i-^i- Then {Si,...,S(i) is a Markov chain 
with 

¥{Sk+i = x\Sk = y) = K-y.p,^,/Y.U,^,pM - y) 

Example 4.2. Let {Ni, . . . ,N(i) ~ iin.{mi....,mi) be a multivariate hyper- 
geometrically distributed random variable, i.e. ¥{Ni = ki,...,Nd = k^) = 

forki G {0,...,m,},...,k, G {0,...,m4 w^th 

ki + . . . + kd = n, and := Ni. Then {Si, . . . , So) is a Markov chain 
with 

F(Sk+i = x\Sk = v) = „^ ^ (x — y) 

5. Definitions and notations for accuracy analyses of algorithms. 

In this section we define terms we need to precisely describe the behaviour 
and the accuracy of numerical algorithms. For M C ]0, oo[ let —M :— {—x : 
X e M} and ±M := M U (-M). 

The lEEE-Double-Precision-Number-System is the set 



IEEE-Double := ±FU±GU {0, -oo, oo} 
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with F := {m • 2^ : m e {2^^ . . . ,2^^ - l},e e {-1074, . . . , 971}} and G := 
{k ■ 2^10^^ : A; e {1, . . . , 2^2 _ ^||^ compare [2]. The values in the defi- 

nition of G and the values (m — 2^^)/2^^ in the definition of F are called man- 
tissas of the considered lEEE-Double-Numbers. We consider the calculation 
of probabilities on computation systems that use lEEE-Double-Precision- 
Numbers. Hence, every computable probability lies in the set IEEE-Double fl 
[0, 1] = G U {m ■ 2^ : m e {2^^ . . . ,2^^ - l},e e {-1074, . . . , -53}} U {0, 1}, 
the minimal computable probability which is greater than zero is min{a; G 
IEEE-Double : x > 0} = 2~^°^^ ^ 5 ■ 10"^^"^ and the maximal computable 
probability which is less than one is max{2; G IEEE-Double : x < 1} = 

1 _ 2-53 ^Q-16_ 

We fix an object not belonging to the set IEEE-Double, call it NaN for 
"Not a Number", and define the four operations 

+,+,-,:: IEEE-Double ^ IEEE-Double U {NaN} 

For x,y e IEEE-Double and o g {+, ■}: 

xoy := min{z G IEEE-Double : z > x oy} 
xoy := max{2; G IEEE-Double : z < x o y} 

except for the following cases: If x = and y G {—00,00} or y = and 
X G {—00, 00} then x~y := x^y := NaN. If x = —00 and y = 00 or y = —00 
and X = 00 then x+y := x+y := NaN. Note that the associative law does 
not hold for these four operations. For example let a = — 1,6 = 1, c = 2~^^, 
then we have a+{b±c) = ^ 2~^^ = {a+b)+c. 

For the calculation of error bounds for the Algorithm A derived in Sec- 
tion |2l we use the following simple fact: 

Lemma 5.1. Let o g {+, ■}, x, y G ]0, oo[ and bi, 62, ci, C2 G IEEE-Double 

with &i < X < Ci and b2 < y < C2- Then 

61062 < X o y < C1OC2 



For a quantitative analysis of the accuracy of computed probabilities we 
need to consider absolute and relative errors. For p,p G [0, 1] we define the 
absolute error 

eabs(p,p) := \p-p\ 

and the relative error 



erei{p,p) := max 



/ eabs(p,p) eabs(l -P, 1 -p) ] ^ \P~P\ 
\ p ^ 1 — p J min(p, 1 — p) 
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in the approximation of p by p, with ^ := and | := oo for x > 0. For 
a, 6 G [0, 1] with a < b and p G [a, b] we further define the absolute error 

eabs([a, b],p) max eabs(p,p) = max{6 -p,p-a} 

pe[a,b] 

and the relative error 

Crei ( [a, 6] , p) : = max e^ei (p, p) 

pe[a,&] 

in the approximation of a probabihty which is known to lie in [a, b] by p. We 
get simple formulas for erei([a, b],p) in the following two cases. If a, 6 G [0,1/2] 
or a, 6 G [1/2, 1] we have 

erei([a,6],p) = max{erei(a,p),erei(6,p)} 

Hence, if a, 6 G ]0, 1/2] we have 

erei([a,o],P) = max{ , — — } 

a b 

and if a, 6 G [1/2, 1[ we have 

/r ,1 cP-a b-p 

erei(a,o,p) =max{- ,- -} 

1 — a 1 — 

For accuracy measurements in interval calculations we use the following mini- 
max errors: 

Definition 5.1. For a,b e [0,1] with a < b we define the absolute 
error 

eabs([a,&]) := min eabs([a, &],p) = eabs([a,6], = 

and the relative error 

erei([a, 6]) := min erei([a, 6],p) 

pe[a,6] 

in the approximation of a probability by the interval [a, b]. 
Easy calculations yield the following formulas: 
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Theorem 5.1. If a,b e [0, 1/2] we have 

2ab 



Vp G [a,b] : erc\{[a,b],p) < e,e\{[a,b], 



a + b 
Hence 

b — a 



erei([a,6]) 

u -r u, 

If a,b E [1/2, 1] we have 

. , a + b — 2ab 

\fp e [a,b\ : erei([a,6J,j9) < erei([a,6J, , 

i/ence 

7) - « 

erei([a,&]) 



2-a-b 



Note that the absolute error eabs(['3,&]) and the relative error erei([a, &]) 
need not be reached simultaneously by one of the approximators. It need not 
be reached at all, as the following example illustrates. 

Example 5.1. In TableUlwe listed the errors eabs([a-, and eT.ei{[a,b],p) 
for [a, b] = [0.02, 0.03] and different approximators p . We see that eabs(['^5 b]) = 
0.005 an c? erei( [a, 6]) = 1/5. If we take the upper houndp = b as approximator 
for the unknown probability p , neither eahsiidy^yp) = eahs{[a,b]) is reached, 
nor erei{[a,b],p) = erei{[a,b]) . 

P eabs([a, e,ci{[a,b],p) 



2a6/(a + 6) = 0.024 0.006 1/5 

(a + 6)/2 = 0.025 0.005 1/4 

a 0.01 1/3 

b 0.01 1/2 

Table 1 



//, for example, the unknown probability is p = (3/10)^ = 0.027, then the 
errors are as listed in Table [B 



We study the maximal accuracy reachable in double-precision probability 
calculations: 
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p 


eahs{p,p) 


eroi(p,p) 


0.024 


0.003 


3/27 


0.025 


0.002 


2/27 


a 


0.007 


7/27 


b 


0.003 


3/27 




Table 2 





Definition 5.2. The maximal accuracy in a double-precision calcula- 
tion of a probability p G [0, 1] is Creiip) '■= erei(/(p)) where I{p) := [max{x G 
IEEE-Double : x < p},mm{x G IEEE-Double : x > p}] is the minimal 
interval containing p, whose endpoints are lEEE-Double-Precision-Numbers. 

Easy calculation yields 

oo < p < 2-1°^^ or 1 - 2-53 < p < 1 
^ m G {1,...,252 - 1}, 

p G 2^^°''"^ ■ ]m, m + 1[ or j9 G 1 — 2-^^ ■ ]m, m + 1[ 



erei(p) 



2m+l 



2^ e {252, . . . ^ - 1}, e G |_io22, . . . , -2}, 

16-52 



p G 2"" 52 . ]^^^ _^ i[ 

p G IEEE - Double n [0, 1] 

From the last formula it follows that 

• in IEEE-Double floatingpoint arithmetic we are able to approximate 
probabilities in [2~^'^'^^, 1 — 2-53j yjQ^ with finite relative error ei.ei(p). 

• for p G [2^^^'^'^, 1/2], the maximal accuracy satisfies erei(p) < 1/(253 _|_ 
1) ^ 1.11 ■ 10-1^ 

6. R Implementation of Markov increment scan algorithms in in- 
terval arithmetic. R is an open source software for statistical computa- 
tions. We extended R by a C-function that, as per C-Standard and lEEE- 
754-Standard [2], allows the operations on IEEE-Double- Numbers which we 
defined in the previous section. We wrote an R-program that implements 
the Algorithm A from Section [2] and uses the principle stated in Lemma 
15.11 to compute bounds for rectangle scan probabilities for Markov incre- 
ments. We implemented the multinomial and multivariate hypergeometric 
transition probabilities, as described in Section |H In a last step the resulting 
R-implementation of Algorithm A sets the returned value to 1, if the original 
return value is greater than 1. 

6.1. Examples. For ~ M„_p with n = 500, d = 365, p = {1 / d, . . . , 1 / d) 
and k E {4, . . . , 32} we computed an upper bound p and a lower bound p for 
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the probability ¥(maxf~^ {Ni + A^j+i + < k) by an R-implementation 

of the Algorithm A from Section |2j In Table |3] we tabulate the computed 
bounds p, p and analyze their accuracy. Numbers written in typewriter font 
are hexadecimal. The coloumn titled "approx" gives the known decimal digits 
of a value of the "probability representation number system" T, that lies 
nearest to the exact value. The probability representation number system T 
consists of all numbers with 7 decimal digits without leading zeros or nines. 
We use the notation .0"^ as an abbreviation for a decimal point followed by x 
zeros, analogously .9^. The symbol ? appearing in a number means that the 
following digits are not exactly known. 

The value Cabs resp. Crei is the minimal upper bound for eabs([P!p]) resp. 
erei([p,p]) which has the form c- lO'^ where c has 3 significant digits and A; G Z. 

Thus, the line with k = 15 means that the probability P(max^~;^^(A^j + 
Ni^i + < 15) lies in the interval [p,p] with 

p = I.fef956911fe58 ■ 2"^ 

= (1 + 15 ■ 16"^ + 14 ■ 16^2 + . . . + 8 ■ 16"^^) ■ 

= 0.99799604913273309847454584087245166301727294921875 

p = 1. f ef 95690c7eda ■ 2~^ 

= (l + 15-16-i + ... + 10-16"^3)-2-^ 

= 0.9979960490927297644958571254392154514789581298828125 

with all equalities exact. The minimal upper bound for eabs(b)P]) which has 
the form c ■ lO'^ where c has 3 significant digits and A; G Z is 2.01 ■ 10^^^ and 
the minimal upper bound for eabs([P;P]) which has this form is 9.99 ■ 10~^. 
A value of the number system T which is nearest to the exact probability is 
0.9979961. As the numbers of the system T in the interval [0.001, 0.9989999] 
differ by 10^'', just knowing the approximate value we can infer that the 
absolute error in this approximation is less than 10"''. 

6.2. Remarks on numerical computations of multinomial probabilities. 

6.2.1. Relative error of complement probabilities. In the preceding sec- 
tion we computed the distribution function of a multinomial scan statistic. 
For several applications, e.g. multinomial scan test, we want to compute the 
upper distribution function instead. If we compute its values from the com- 
plements P(maxt 1 Ni + Ni+i + Ni+2 >k) = l- P(maxt ^ Ni + A^i+i + Ni+2 < 
k — 1) in exact arithmetic, for example by using a suitable software, there 
is no increse of error. If we do automatic computation of the complement in 
lEEE-Double-Precision-Number-System, then the error increases for small 
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k 


P,P 


^abs 




approx 


4 














5 


I.lc5dflelalf83 ■ 2''-''^ 
I.lc5dflel71043 ■ 2-^'^* 


5.82 ■ W-'^^ 


2.01 ■ IQ-ii 


.0^^28993 


6 


I.b826f22f 10057 ■ 2-"' 
I.b826f22ec43c3 ■ 2'^''^ 


2.34 ■ lO^^i 


2.01 ■ 10-11 


. 01^11651 


7 


I.b71c492587c97 ■ 2"^' 
I.b71c49253c2df ■ 2'^'' 


2.57 ■ 10-1^ 


2.01 ■ 10-11 


.0''12780 


8 


1.98b8351d76fbd ■ 2"^^ 
1.98b8351d309cf ■ 2"" 


1.57 ■ 10^1'' 


2.01 ■ 10-11 


.0^77957 


9 


1.0f0230ce6f8al • 2"" 
1.0f0230ce40el5 • 2"* 


1.33 ■ 10^12 


2.01 ■ 10-11 


.0661642 


10 


1.826e2adb7befd • 2"^ 
1 . 826e2adb39686 • 2"^ 


7.57 ■ 10-12 


2.01 ■ 10-11 


.3773734 


11 


1.7131cf887a229 • 2^1 
1.7131cf883a935 • 2"^ 


1.45 ■ 10-11 


5.19 ■ 10-11 


.7210832 


12 


1 . ce576094ddb84 • 2"^ 
I.ce5760948elf6 • 2"^ 


1.81 ■ 10-11 


1.87 ■ 10-1'^ 


.9030104 


13 


l.f 1162301d80ec • 2"^ 
l.f 1162301827ae • 2"^ 


1.95 ■ 10-11 


6.69 ■ 10-1'^ 


.9708720 


14 


I.fbef9498b0df9 • 2"^ 
I.fbef9498596d7 • 2"! 


1.99 ■ 10-11 


2.51 ■ 10-9 


.9920622 


15 


I.fef956911fe58 • 2"^ 
I.fef95690c7eda • 2"i 


2.01 ■ 10-11 


9.99 • 10-9 


.9979961 


16 


I.ffclfbbfd6e58 ■ 2"^ 
l.ffclfbbf7ecbl ■ 2"^ 


2.01 ■ 10-11 


4.24 ■ 10-* 


.9^52685 


17 


I.fff23b0d23a3c ■ 2"^ 
I.fff23b0ccb810 ■ 2'^ 


2.01 ■ 10-11 


1.91 ■ lO-'^ 


.9^89495 


18 


I.fffdld22cb527 • 2"^ 
l.fffdld22732da ■ 2"^ 


2.01 ■ 10-11 


9.11 ■ lO-'^ 


.9*77980 


19 


I.ffff6d5024936 ■ 2"^ 
I.ffff6d4fcc6e4 ■ 2"! 


2.01 ■ 10-11 


4.59 • 10-^ 


.9^56284 


20 


I.ffffe4570f39a • 2'^ 
I.ffffe456b7146 • 2"^ 


2.01 ■ 10-11 


2.44 ■ lO-'"^ 


.9«'17567 


21 


I.fffffb08bdl3c ■ 2-^ 
I.fffffb0864ee9 ■ 2"! 


2.01 ■ 10-11 


1.36 ■ 10-1 


.9«'8520? 


22 


I.ffffff264f47d • 2-^ 
I.ffffff25f7228 ■ 2-1 


2.01 ■ 10-11 


7.91 ■ 10-1 


.9'^74? 


23 


l.ffffffdc79315 ■ 2-1 
I.ffffffdc210c0 ■ 2-1 


2.01 ■ 10-11 


4.83 • 10-3 


.9*6? 


24 


l.fffffffa913ba ■ 2'^ 
l.fffffffa39167 ■ 2-1 


2.01 ■ 10-11 


3.08 ■ 10-2 


.9"? 


25 


I.ffffffff53a50 • 2-1 
l.fffffffefb7fe • 2^1 


2.01 ■ 10-11 


2.04 ■ 10-1 


.99? 


26 


1 

I.ffffffffb44b7 ■ 2-1 


1-p 


oo 


.91"? 


27 


1 

l.ffffffffcf373 ■ 2^1 


1-p 


oo 


.91"? 


28 


1 

I.ffffffffd2fd3 • 2^1 


1-p 


oo 


.91"? 


29 


1 

l.ffffffffd37fa ■ 2-1 


1-p 


oo 


.91"? 


30 


1 

l.ffffffffd3908 • 2^1 


1-p 


oo 


.91"? 


31 


1 

l.ffffffffd392a • 2^1 


1-p 


oo 


.91"? 


32 


1 

l.ffffffffd392a • 2^1 


1-p 


oo 


.91"? 



Upper and lower bounds p,p for ¥{Tiiaxf^^ Ni + A^^+i + iVi+2 < k) with N ^ M„^p, 
n = 500, d = 365, p — (l/d, . . . , l/d) and k £ {4, . . . , 32}. For details, see Suhsection \6.1\ 
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probabilities. Then, we are not able to approximate probabilities less then 
10~^^ with a finite relative error. Compare Table |H Here, the relative er- 
ror increases for small probabilities as well as for big probabilities. Small 
complements of probabilities are lost. In general, one should try to avoid de- 
veloping algorithms that complement the computed probability at the end. 
An algorithm that computes the complement is not equivalent to a direct 
one. 

The maximal accuracy with respect to complementation of probabilities is 
defined as 



ereiAp) '■= max(erei(p),erei(l - p)) 



Easy calculation yields 

{oo < p < 2"^^ or 1 - 2~^^ < p < 1 
^ mG{l,...,2^^-l}, 
p e 2"^^ ■]m,m + l[oT p e 1 - 2~^^ ■ ]m, m + 1[ 
pe IEEE - Double n [0, 1] 

6.2.2. Computation Time and Space. Besides the accuracy of the algo- 
rithm, there are two other problems that matter: Time and space needed to 
compute the probability. 

The implementation of the multinomial scan algorithm we made needs to 
store 2 * ("^^) Double-Precision-Numbers. Each Double-Precision-Number 
needs 8 Bytes. For example, for the scan width £ = 3, on a computer with 
16 GByte memory, we were able to compute Scan-Probabilities for up to ap- 
proximately n = 1700 in double-precision. Using Single-Precision-Numbers, 
which take only 4 Bytes, we could compute up to n = 2150, but the accuracy 
is worse than in double-precision computations, as Table O demonstrates. 
The lEEE-Single-Precision-Number-System is the set 

IEEE-Single := ±F U ±G U {0, -oo, oo} 

with F := {m*2^ : m e {2^3, . . . , 2^^ - 1}, e G {-149, . . . , 104}} and G : = 
{k * 2~^^^ : A; G {1, . . . , 2^'^ — 1}}, compare [2]. In the third coloumn of Ta- 
ble [5] we listed the first digits of the computed bounds p,p in decimal format. 

The time that it takes to compute a rectangle scan probability for a multi- 
nomially distributed random vector in single precision does not differ much 
from the time it takes in double-precision, examples are listed in Table Ei 



6.3. Binomial Probabilities. We use the following algorithm to compute 
the multinomial transition probabilities, that are binomial. 



RECTANGLE SCAN PROBABILITIES FOR MARKOV INCREMENTS 
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k 


P,P 


^abs 




approx 


5 


1 








1 


6 


1 

l.fffffffffffff ■ 2^1 


l-£ 


oo 


.9"? 


7 


1 

l.fffffffffffff ■ 


1-p 


oo 


.9"? 


8 


l.ffffff9238edc • 2"^ 
l.ffffff9238edb ■ 2'^ 


5.55 • 10"" 


4.34 • 10-8 


.9'^87220 


9 


I.ff99dlf2b8b3e • 2"^ 
I.ff99dlf2b8a24 • 2"^ 


1.57- 10-1"* 


2.01 ■ IQ-ii 


.9322042 


10 


I.delfb9e637e3e • 2"^ 
I.delfb9e6320eb • 2"! 


1.33 ■ 10-12 


2.01 ■ IQ-ii 


.9338358 


11 


1.3ec8ea92634bd • 2"^ 

1.3ec8ea9242081 ■ 2~l 


7.57-10-12 


2.01 ■ 10-11 


.6226266 


12 


I.ld9c60ef8ad96 • 2"^ 
I.ld9c60ef0bbae ■ 2"^ 


1.45 ■ 10-11 


5.19 ■ 10-11 


.2789168 


13 


1.8d44fb5b8f050 • 2"" 
1.8d44fb59123e0 ■ 2-" 


1.81 • IQ-ii 


1.87- 10-1° 


.0969896 


14 


I.dd3b9f Cfb0a40 • 2"" 
I.dd3b9f c4fe280 ■ 2"" 


1.95 • IQ-ii 


6.69 ■ 10-1° 


.0291280 


15 


1 . 041ad9e9a4a40 ■ 2" ' 
1.041ad9d3c81c0 ■ 2"^^ 


1.99 • IQ-ii 


2.51 • 10-8 


.0079377 


16 


1.06a96f 3812600 • 2-"' 
1.06a96ee01a800 • 2 '' 


2.01 • 10 11 


<).!)!) ■ 10 " 


,0020040 


17 


I.f0220409a7800 • 2"^^ 
I.f0220148d4000 • 2'^"^ 


2.01 ■ IQ-ll 


4.24 • 10-8 


.0^47315 


18 


I.b89e668fe0000 • 2-^" 
I.b89e5b8b88000 ■ 2"^"' 


2.01 • IQ-ll 


1.91 • 10-''' 


.0^10505 


19 


1.716ec66930000 • 2"^" 
1.716e9a56c8000 ■ 2^1^ 


2.01 • IQ-ll 


9.11 • 10-''' 


.0^22020 


20 


1.2560672380000 • 2"^" 
1.255fb6d940000 ■ 2-^** 


2.01 ■ IQ-" 


4.59 • 10-8 


.0^4371? 


21 


I.ba948eba00000 • 2"^^ 
l.baSfOceeOOOOO ■ 2-21 


2.01 • 10-11 


2.44 • 10-5 


.0^824? 


22 


1.3de6c45c00000 ■ 2 ^ ^ 
1.3dd0bbl000000 -2-23 


2.01 ■ 10-11 


1.36 ■ 10-"* 


.0'*14? 


23 


I.b411bb0000000 • 2''^'^ 
l.b361706000000 • 2-^6 


2.01 ■ 10-11 


7.91 ■ 10-* 


.0''253? 


24 


1 . lef 7a00000000 • 2-'-*** 
l.lc36758000000 ■ 2^28 


2.01 ■ 10-11 


4.83 ■ 10-2 


.0^41? 


25 


1.71ba640000000 • 2-'^^ 
1.5bbl 180000000 ■ 2-''''^ 


2.01 ■ 10-11 


3.08 ■ 10-2 


.086? 


26 


1.0480200000000 • 2-^^ 
1 . 58b6000000000 • 2-3* 


2.01 ■ 10-11 


2.04 • 10-1 


.0''? 



Table 4 

Upper and lower hounds p,p for P(max^~^^ Ni + A^i+i + > k) with N ^^M. 
n = 500, d = 365, p={l/d,..., 1/d) andkG{5,..., 26}. 
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k 






^abs 




approx 


4 

















5 


1.974c00 • 2-"^^ 



.0"'"3652... 




1.83 • lO-^i 


1.04 • 10-2 


.O'lO? 


6 


I.bcc5a4 • 2"'" 
l.b39300 ■ 2-S^ 


.0i''1177... 
.0i''ll52... 


1.22 • 10-22 


1.04 ■ 10-2 


.Oi^ll? 


7 


l.bbb862 ■ 2"2v 
I.b28b40 -2^27 


.0''12913... 
.0^12646... 


1.34 • 10-1° 


1.04 • 10-2 


.o^'m 


8 


1.9d02a2 ■ 2"^^ 
1.947834 • 2-" 


.0^78775... 
.0377146... 


8.15 ■ IQ-f' 


1.04 ■ 10-2 


.037? 


9 


l.llda84 • 2-" 
l.OcSOdO • 2-* 


.0668587... 
.0654762... 


6.91 ■ 10-4 


1.04 ■ 10-2 


.06? 


10 


1.867cac • 2-^ 
1.7e699a • 2"^ 


.3813349... 

.3734497... 


3.94 ■ 10--^ 


1.04 • 10-2 


.3? 


11 


1.7511f c • 2-^ 
1.6d5b2a • 2"^ 


.7286528... 
.7135861... 


7.53 ■ 10-'^ 


2.7- 10-2 


.7? 


12 


I.d331e6 • 2-^ 
l.c988cc ■ 2-1 


.9124900... 
.8936218... 


9.43 • 10-3 


9.7- 10-2 


.? 


13 


I.f64e04 • 2-^ 
l.ebebie • 2-^ 


.9810639... 
.9607779... 


1.01 • 10-2 


3.49 ■ 10-1 


.9? 


14 


1 

I.f6a7a6 • 2-^ 


1 

.9817478... 


l-p 


oo 


.9? 


15 


1 

I.f9a956 • 2-1 


1 

.9876200... 


l-p 


oo 


.9? 


16 


1 

I.fa6fe6 • 2-1 


1 

.9891349... 


l-p 


oo 


.9? 


17 


1 

I.fa9fa0 • 2-1 


1 

.9894990... 


l-p 


oo 


.9? 


18 


1 

l.faaa68 ■ 2-1 


1 

.9895813... 


l-p 


oo 


.9? 


19 


1 

l.faacbe • 2-1 


1 

.9895989... 


l-p 


oo 


.9? 


20 


1 

l.faad2c ■ 2-i 


1 

.9896024... 


l-p 


oo 


.9? 


21 


1 

l.faadSc • 2-1 


1 

.9896029... 


l-p 


oo 


.9? 


22 


1 

l.faad40 • 2-1 


1 

.9896030... 


l-p 


oo 


.9? 


23 


1 

l.faad44 • 2-1 


1 

.9896031... 


l-p 


oo 


.9? 


24 


1 

l.faad46 • 2-1 


1 

.9896032... 


l-p 


oo 


.9? 


25 


1 

l.faad46 • 2-1 


1 

.9896032... 


l-p 


oo 


.9? 



Table 5 

Upper and lower hounds p,p for P(niax^J^^ Ni + iV^+i + < k) with N ~ M„^p, 
= 500, d = 365, p = {1/d, . . . ,1/d) and fc € {4, . . . , 25}, computed in single-precision. 
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n 


d 


k 


p (double) 


time (double) 


p (single) 


time (single) 


100 


365 


6 


0.9934578 


1 s 


0.9914927 




1 s 


500 


365 


13 


0.9708720 


1 min 27 s 


0.9607779 


1 


min 25 s 


1000 


365 


20 


0.9604324 


49 min 39 s 


0.9405573 


10 


min 19 s 


1500 


365 


27 


0.9739303 


3 h 22 min 26 s 


0.9438554 


2 h 44 


min 00 s 


1700 


365 


29 


0.9610315 


5 h 47 min 01 s 


0.9274842 


4 h 58 


min 21 s 


1750 


365 


31 


X 


no computation possible 


0.9516879 


6 h 36 


min 48 s 


2150 


365 


37 


X 


no computation possible 


0.9507257 


12 h 47 


min 22 s 



Table 6 

Computation time for a lower bound p for the scan probability 



P(max^^;^^ Ni + Ni^i + -/Vi-|_2 < k) for a random vector (iVi, . . . , Nd) with the multinomial 
distribution Mn,p with p = (1/rf, ■ ■ ■ T^/d), in single precision and in double precision. 
Details are described in Subsubsection \E. 2. 2\ 

double bnp (unsigned int k, unsigned int n, double p, double q){ 
if (2*k>n) return(bnp(n-k,n,q,p)) ; 
double f=1.0; 

unsigned int jO=0, jl=0, j2=0; 

while ( (jO<k) I (jl<k)| (j2<n-k) ) 

{ 

if( (jO<k) && (f<l) ) {jO++; f*= (double) (n-k+jO)/ (double) jO;} 
else { if(jl<k) {jl++; f*= p;} else {j2++; f*= q;} } 

} 

return f ; 
} 

For upper bounds p and q and "rounding up" mode the algorithm calculates 
an upper bound for the exact binomial probability. For lower bounds p and 
q and "rounding down" mode the algorithm calculates a lower bound for the 
exact binomial probability. 

6.4. Hypergeometric Probabilities. We use the following algorithm to com- 
pute the multivariate hypergeometric transition probabilities, that are uni- 
variate hypergeometric. Table [7] contains the distribution function of the 
random variable max'fl^{Ni + A^j+i + with ~ H„^m with n = 500, 

rf = 365 andm= (10, ...,10). 

double hyp (int n, int r, int b, int k){ 

double f=1.0; 

int jO=0, jl=0, j2=0; 

while ( (jO<k)| (jl<n-k) I (j2<n) ){ 

if(f<l && ( (jO<k) I (jl<n-k)) ){ 

if (jO<k) { f*=(double) (r-jO)/(jO+l) ; jO++;} 

else {if (jl<n-k) { f *=(double) (b-jl)/(jl+l) ; jl++;} 
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else if (j2<n) -[f *=(double) (r+b-j2)/(j2+l) ; j2++;}} 
> 

else if (j2<n) i f*=(double) (j2+l)/(r+b-j2) ; j2++;} 
} 

return f ; 
} 

In the "rounding up" mode the algorithm calculates an upper bound for the 
exact hypergeometric probability. In the "rounding down" mode the algo- 
rithm calculates a lower bound for the exact hypergeometric probability. 
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k 


[P,P] 


^abs 


^rel 


approx 


4 














5 


1.94a78cce6bf78 • 2"^"" 
1.94a78cce088a0 ■ 2-"'" 


3.09 ■ 10-5^ 


2.86 ■ 10^11 


.0*^10815 


6 


1.0acc3dae78827 ■ 2"-'-' 
1.0acc3dae36d0e ■ 2^^'^' 


8.29 • IQ-^® 


2.87 • IQ-ii 


.01^^28926 


7 


1.591d6928456d6 ■ 2"^" 
1.591d6927f05d0 ■ 2^20 


3.69 • IQ-i''' 


2.87 • IQ-ii 


.0^12856 


8 


1.40ac4ad3593a9 -2-' 
1.40ac4ad30a26f -2"^ 


2.81 • IQ-i^ 


2.87 • IQ-ii 


.0097862 


9 


I.df885f4b6ceae • 2"-' 
I.df885f4af6a55 • 2 


1,1)1 • 10 " 


2.87 • in ^-^ 


.2341 168 


10 


1.546bd869a7f5e • 2"^ 
1.546bd86953fe9 • 2"! 


2.60 • 10-11 


5.70 • IQ-ii 


.6648853 


11 


I.ceclebd5b5793 • 2~^ 
l.ceclebd543545 • 2"^ 


2.81 • 10-11 


2.70 • IQ-i" 


.9038233 


12 


I.f4e8088a29393 • 2"^ 
I.f4e80889adab5 • 2^^ 


2.86 • 10-11 


1.30 • 10-® 


.9783328 


13 


I.fde26f4234a4c • 2"^ 
I.fde26f41b6dfc • 2^^ 


2.87 • 10-11 


6.92 • 10-® 


.9958682 


14 


I.ffa6780ca228e • 2"^ 
I.ffa6780c23f48 ■ 2^1 


2.87 • 10-11 


4.20 • 10-* 


.9^31693 


15 


I.fff314a41d498 • 2"^ 
I.fff314a39f023 • 2"! 


2.87 ■ 10^11 


2.91 • IQ-'^ 


.9*01433 


16 


I.fffe5ec7c001c • 2"^ 
I.fffe5ec741b7c • 2"! 


2.87 ■ IQ-ii 


2.31 ■ lO"" 


.9*87566 


17 


I.ffffd2049693f • 2"^ 

l.ffffd20418497 • 2"^ 


2.87 ■ 10^11 


2.10 ■ ID"-'' 


.9^8629? 


18 


l.fffffb9535338 • 2"^ 
I.fffffb94b6e8e • 2"^ 


2.87 ■ 10^11 


2.18 ■ 10^'' 


.9^868? 


19 


I.ffffffaldc0a3 • 2"^ 
l.ffffffalSdbfb • 2-^ 


2.87 • 10-11 


2.61 • 10-3 


.9''8? 


20 


l.fffffff9717bl • 2-^ 
I.fffffff8f330a • 2^1 


2.87- 10-11 


3.63 • 10-2 


.9»? 


21 


l.ffffffffd3bf 1 • 2-^ 
l.ffffffff 55749 • 2-1 


2.87- 10-11 


5.88 ■ 10-1 


.910? 


22 


1 

l.ffffffffbb782 • 2-1 


1-p 


oo 


.910? 


23 


1 

I.ffffffffc0de3 • 2-1 


1-p 


oo 


.910? 


24 


1 

l.ffffffff cllb4 • 2-1 


1-p 


oo 


.910? 


25 


1 

l.ffffffff clld9 • 2^1 


1-p 


oo 


.910? 


26 


1 

l.ffffffff clld9 • 2-1 


1-p 


oo 


.910? 



Table 7 

Upper and lower hounds p,p for F{-niaxf~^ Ni + Ni^i + iVi+2 < k) with A'' ~ H, 
n = 500, d = 365, m = (10, . . . , 10) and k G {4, . . . , 26}. 



